Controlled Markovian dynamics of graphs: 

unbiased generation of random graphs with 

prescribed topological properties 
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Abstract — We analyze the properties of degree-preserving 
Markov chains based on elementary edge switchings in undi- 
rected and directed graphs. We give exact yet simple formulas 
for the mobility of a graph (the number of possible moves) 
in terms of its adjacency matrix. This formula allows us to 
define acceptance probabilities for edge switchings, such that the 
Markov chains become controlled Glauber-type detailed balance 
processes, designed to evolve to any required invariant measure 
(representing the asymptotic frequencies with which the allowed 
graphs are visited during the process). 

I. Introduction 

Sampling uniformly the space of graphs with prescribed 
macroscopic properties has become a prominent problem in 
many application areas where tailored graph ensembles are 
used as proxies or null models for real complex networks. 
Quantities measured in real networks are often compared 
with the values these quantities take in their randomised 
counterparts. Ensembles of randomised networks allow one to 
put error bars on these values, and therefore to identify which 
topological features of real networks deviate significantly from 
the null model. Such features are likely to reflect e.g. design 
principles or evolutionary history of the network. Ensembles 
of randomised graphs are to be generated numerically, and 
each graph realisation should be produced with a probability 
proportional to a prescribed statistical weight (often taken to 
be uniform) of the graph, for the analysis to be unbiased. 
At the level where the constraints of the randomized graphs 
involve only the simplest quantities, i.e. the degree sequence, 
even uniform generation of such graphs is known to be a non- 
trivial problem fl]-||4). One classical algorithm for generating 
random networks with prescribed degrees [5 j assigns to each 
node a number of 'edge stubs' equal to its desired degree, 
and joins iteratively pairs of randomly picked stubs to form a 
link. A drawback of this algorithm is the need for rejection of 
forbidden graphs (those with multiple edges or self- loops), 
which can lead to biased sampling [6|. A second popular 
method for generating random graphs with a given degree 
sequence is 'edge swapping', which involves successions of 
ergodic graph randomising moves that leave all degrees in- 



variant Q, JS). However, naive accept-all edge swapping will 
again cause sampling biases. The reason is that the number of 
edge swaps that can be executed is not a constant, it depends 
on the graph c at hand; graphs which allow for many moves 
will be generated more often. Any bias in the sampling of 
graphs invalidates their use as null models, so one is forced to 
mistrust all papers in which observations in real graphs have 
been tested against null models generated either via the 'stubs' 
method or via randomisation by 'accept-all edge swapping'. 
This situation can only be remedied via a systematic study of 
stochastic Markovian graph dynamics, which is the topic of 
this paper. We determine the appropriate adjustment to the 
probability of accepting a randomly chosen proposed edge 
swap for the process to visit each graph configuration with 
the same probability. This can be done for nondirected and 
directed graphs. We will also show that our method can be 
used to generate graphs with prescribed degree correlations. 

II. Generating random graphs via Markov chains 

A general and exact method for generating graphs from the 
set G[k] — {c e G\ k(c) = k] randomly (where k denotes the 
degree sequence, or the joint in- and out-degree sequence of 
the graph), with specified probabilities p(c) = Z l exp[-H(c)] 
was developed in |9|. It has the form of a Markov chain: 

VceG[k]: p t+l (c) = £ W(c\c')p t (c') (1) 

C'eC[k] 

Here p t (c) is the probability of observing graph c at time t in 
the process, and W(c\c') is the one-step transition probability 
from graph c' to c. For any set <t> of ergodic reversible 
elementary moves F : G[k] — > G[k] we can choose transition 
probabilities of the form 



W(c\c') = J]q(F\c') 



6 c ,fc'A(Fc'\c') + S c ,c'[l-A(Fc'\c')] 

(2) 



The interpretation is as follows. At each step a candidate move 
F € <D is drawn with probability q(F\c'), where c' denotes the 
current graph. This move is accepted (and the move c' — > 



c - Fc' executed) with probability A(Fc'\c') 6 [0, 1], which 
depends on the current graph c' and the proposed new graph 
Fc'. If the move is rejected, which happens with probability 
l-A(Fc'\c'), the system stays in c'. We may always exclude 
from O the identity operation. One can prove that the process 
(Q3 will converge towards the equilibrium measure pco(c) - 
Z _1 exp[-i/(c)] upon making in (f2]l the choices 



q(F\c) 
A(c\c') 



I F (c)/n(c) 



n(c , )e -\[H(C)-H(C')] 



(3) 
(4) 



n(c , )e -{[H(c>-H(C')] + n^imo-HiC')] 

Here /f(c) = 1 if the move F can act on graph c, If(c) - 
otherwise, and n(c) denotes the total number of moves that 
can act on a graph c (the 'mobility' of state c): 



n(c) = ^ I F (c). 



(5) 



FeO 



III. Degree-constrained dynamics of nondirected graphs 



We first apply our results to algorithms that randomise 
undirected graphs, while conserving all degrees, by application 
of edge swaps that act on quadruplets of nodes and their 
mutual links. Such moves were shown to be ergodic, i.e. any 
two graphs with the same degree sequence can be connected 
by a finite number ofsuccessive edge swaps 1171, Il8l. 

Let us define the set Q - {(i,j,kj) e {l,...,N} 4 \ i<j<k< 
() of all ordered node quadruplets. The possible edge swaps to 
act on (i, j, k, () are the following, with thick lines indicating 
existing links and thin lines indicating absent links that will 
be swapped with the existing ones, and where (IV, V, VI) are 
the inverses of (I, II, III): 
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We group the edge swaps into three pairs (I,IV), (II,V), and 
(III, VI), and label all three resulting auto-invertible operations 
for each ordered quadruple (i,j,k,£) with a subscript a. Our 
auto-invertible edge swaps are now written as Fjj k ( ;a , with 
i < j < k < ( and a e {1,2,3}. We define associated indicator 
functions hjki^ic) € {0, 1} that detect whether (1) or not (0) 
the edge swap Fy# ;o can act on state c, so 

hjkt-\{c) = CijC ke (l-Cic)(l-c kj ) + (1 -Cjj)(\-c kl )Citc k j (6) 
hjkt-iic) = CijC lk (\-ci k )(\-c (] ) + (1 -Cjj)(\-C( k )c ik C(j (7) 
hjm(c) = c ik Cji{\-Ci(){\ -c jk ) + (1 -c ik ){\ -Cji)c it Cj k (8) 



If hjkt;a(c) — 1> this edge swap will operate as follows: 

Fij k (- a (c)gr = 1 ~C qr for (q, f) 6 SijU;a (9) 

Fij k t- a {c) qr - c qr for (q,r)iSfjU;a (10) 

where 

Siju-i = {(i, j), (*,0, ft 0, (*,/)} (11) 

Stjua = (ft /),(**). ft *),(*,/)} (12) 

Stjua = {ft *),U <),(*,*),(/,*)} (13) 

Insertion of these definitions into the recipe (121314b gives 

Iijk€;a(c') 



W(c\c') = 2 Z 

i<j<k<( ff<3 



n(c') 



(14) 



3C,F iMa ,C 



, s -\[E(F ljkc , a C')-E{C')} + g c c ,^lE(F ukr , a C')-E(C')] 



Q'\[E{F ijk[ ,„C')-E{C')] + Q \[E{F iju , a C')-E{C')] 

with E(c) - H(c)+logn(c). The process ([141 can be described 
as the following algorithm. Given an instantaneous graph c'\ 
(i) pick uniformly at random a quadruplet (;', j, k, €) of sites, (ii) 
if at least one of the three edge swaps c' — > c - Fij k ( ia (c') is 
possible, select one of these uniformly at random and execute 
it with an acceptance probability 



A(c\c') = [l +e E W«^0-E(co]- 



(15) 



then return to (i). For this Markov chain recipe to be practical 
we finally need a formula for the mobility n(c) of a graph. 
This could be calculated [9], givingj (with TrA = 2/^;;) : 



n(c) 



- iK*0 2+ zE*-5E*-5l><A 



+ lTr(c 4 )+iTr(c 3 ) 



(16) 



Naive 'accept-all' edge swapping, where A(c\c') - 1, corre- 
sponds to choosing E(c) — in ( TBI , and would give the biased 
graph sampling probabilities poo(c) - n(c)/ Y>c' n ( c ') upon 
equilibration. The graph mobility acts as an entropic force 
which can only be neglected if (fTol l is dominated by its first 
three terms; it was shown that a sufficient condition for this to 
be true is {k 2 )k max /{k) 2 <s:N. In networks with narrow degree 
sequences this condition holds, and naive edge swapping is 
roughly acceptable. However, one has to be careful with scale- 
free graphs, where (k 2 ) and £ max diverge as A^ —> oo. 

IV Degree-constrained dynamics of directed graphs 

For directed networks the generalisation of the canonical 
edge swap (up to relabelling of nodes) is shown below: 




1 Efficient expressions for the change in mobility following one move have 
also been derived |9], to avoid unnecessary matrix multiplication in the 
computational implementation. 



This move samples the space of all directed graphs with 
prescribed in- and out- degrees ergodically only if self- 
interactions are permitted [10]. We could now proceed with 
our formalism as before, but it is more efficient in the case of 
directed graphs to define the swaps and the indicator function 
in terms of pairs of links rather than quadruples of nodes. 

Let c now be a nonsymmetric connectivity matrix, such 
that c a i, — 1 if a — > b, otherwise c a b - 0. Consider A to be 
the set of links within the network defined by c. We write 
x = (Xi,Xj) € A if and only if c XiX . = 1. The indicator function 
can be written as 



'x,y;n 



1 if x, y E A and (x h yj), (y h Xj) i A 
otherwise 



If ^x,y;n = 1, then the corresponding autoinvertible operation 
Fx,y,n acts on state c as follows: 

F x ,y,n(c) qr = 1 - c qr for q e {*;,)>,-} and r 6 {Xj,yj} (17) 
Fx,y-n(c)q r = c qr otherwise (18) 

When self-interactions are forbidden, a further type of ele- 
mentary move is required to ensure ergodicity [10]. It can be 
visualised as reversing a triangle cycle: 

yj jj 




Xi 



The indicator function for this new move is 



'x,y;A 



1 if x, y, (y h xd e A and Xj = y ; 

x-\y- l ,(xj,yj) £A 
otherwise 



where x = (xj, xi) represents a link that is to be reversed. The 
corresponding autoinvertible operation represented by Fx,y,& 
acts on state c as follows: 



r x,y;A\C)qr 
r x,y;A\C)qr 



1 



- qr 



- qr 



for(q,r)eS Xl , Xi , yi (19) 
for (q,r) t S Xi , Xj , y} (20) 



in which S a b c - {(a,b),(b,c),(c,a),(b,a),(c,b),(a,c)} is the 
set of pairs of vertices of the triangle. 

We now combine edge swaps and triangular reversions into 
(0. This requires us to express 2fe<D l(F\c') from the point of 
view of sampling links x and y. This is in fact straightforward, 
as each pair of links can either trigger the square indicator 
function (if they define an edge swap), a triangle indicator 
function (if they define a triangle reversion), or neither (when 
the indicator function returns zero). It only remains to observe 
the symmetries of the problem: cycling through each pair of 
bonds we will come across each unique edge swap twice, and 
across each unique triangle swap three times IfTTI . Hence 



2>ev)=2 



2-'x,y;n t 3 'x,y;A 



x.yeA 



n(c') 



which expresses sampling over all possible moves in terms 
of (the more tangible) sampling over links. With the above 



definitions in hand, the transition probability is given by IfTTI 
W(c\c>) = X ^j^ (21) 



x.yeA 



n(c') 



5c f* r e-J [£(F ^ c,) - £(C ' )] + S c cei lE(F *> c ' yE{C ">h 



e -{[E(F SJ C')-E(C')] + e {[E(F K jC')-E(C')] 

where n(c) - n n (c) + n A (c) is the total number of valid moves 
that can act on the directed graph c, written as the sum of the 
number of edge swaps and the number of triangle reversions. 
The mobility terms can again be calculated, giving IfTTI 

n n (c) = -Tr(cc T cc T ) -^ k°%kj + Tr(cc T c) 

+^y>) 2 + irr(c 2 ) -Vkfk" 1 



n A (c) 



I I 

-Tr(c 3 ) - 7Y(c : c 2 ) + 7Y(c l2 c) - -7Y(c :3 ) 



with (c T )jj = Cji and (c*)y = c/jCj/. Also for directed graphs 
it is possible to calculate efficient expressions for how the 
mobilities change following a single move. 

For the simplest case where H{c) is constant, i.e. the Markov 
chain is to evolve towards a flat measure, we immediately 
observe that the required move acceptance probabilities are 

At i >\ h ^ n n (c) + n A (c) r i 

A(c\c) = 1 + — r-r — (22) 

L n D (c') + nA(c') J 



V. Numerical examples for directed graphs 
Examples illustrating the canonical Markov chains for 
nondirected graphs are given in (9). Here we compare two 
variants of the directed edge rewiring algorithm: the naive 
'accept all moves' version, and the version with the canonical 
mobility corrections (f22l . Consider a network with K + 2 
nodes, of which one has degrees (& m ,£ out ) = (0, K), one 
has (k m ,k° ut ) = (A^ 0), and the remaining K have degrees 
(£ in ,£: out ) = (1,1). Self-interactions are forbidden. Up to 
relabelling, there are two such graph types (see below): 




K=25 


accept all 
process 

avg. n(c) 


correct 
process 

avg. n(c) 


predict 
actual 


58.52 
58.32 


47.92 
47.95 



The left graph type has mobility K(K- 1), and only occurs in 
one network configuration. The right graph type has mobility 
of 2K-3, and multiplicity K(K-l). The table shows the results 
of numerical simulations for K = 25, compared to theoretical 
predictions. We have used the mobility itself as a marker for 
the proportion of time spent in each type of configuration; 
the difference between the two processes is striking, and will 
translate into differences in measurements of any observable 
which differentiates between the two configurations. 



P(k) 




n(*,tf|c ) 

* 




LT(fc, k') (theory) 



n(A:,fc'|c fl „ai) 




52=1 




Fig. 1. Results of Markovian dynamics tailored to target the non-uniform 
measure J23I . Top left: degree distribution of the (randomly generated) initial 
graph c , with N = 4000 and (k) = 5. Top right: II(k,k'\c ) of the initial 
graph. Bottom left: the target relative degree correlations {26) chosen in )23t . 
Bottom right: colour plot of TI(k, &'|eflnal) in the final graph Cfi na |, measured 
after 75,000 accepted moves of the Markov chain U4t . 

VI. Generation of random graphs with prescribed degree 

CORRELATIONS VIA REWIRING ALGORITHMS 

Our approach can be extended to accurately target desired 
degree correlations defined by W{k, k'). This can be achieved 
for undirected graphs by ensuring convergence of the above 
Markov chain to the following non-uniform measures iTTZl . 

(23) 



Pic) 



n 



N p(ki)p{kj) 



<V + (i- 



n p(h: 



with p(k) = N^Zid^, W(k,k') = (N(k))- 1 l, i jc tj 6 kl jjk,j lf , 
<*> = !**/>(*), and 

W(jfc) = ^ W(fe, /f ) = p(k)k/(k) (24) 

In the language of our process, (l23l corresponds to the choice 



ff(c) = -2l°g 



'<] 



Uk) W(kj,kj) 

[ N pikdpikj) 



*«„1 + 



(- 



<jt> W(fc,*/) 



A? p(*i)j>(* 



^K° 



Hence, for the candidate edge swaps c' — > c — Fijkc-, a c' the 
acceptance probability (@]i can be used, where 



H(C)-H(C) 






(25) 



with L a b - Nk /[Tl(k a , kb)k a k D ] - 1 and relative degree correla- 
tions Tl(k,k') = W(fc,*0/W(ifc)W(if). 

We show an example of the Markov chain (fl4l targeting 
(1231 where relative degree correlations are chosen as 

nOfc,Jf) = (k-k'fllPx-fak+p^m-fak' +/3 3 k' 2 ](26) 



(the parameters /?,■ follow from d24l >). An initial graph cq was 
constructed with a non-Poissonian degree distribution and no 
degree correlations, i.e. Yl(k, k'\co) ~ 1. After iterating the 
Markov chain until equilibrium (after 75, 000 accepted moves, 
and after reaching maximal Hamming distance between initial 
and final configuration) degree correlations are seen in very 
good agreement with their target values; (see Figure [T|. 
Extension to directed graphs is achieved by replacing k with 
k = (£ ln ,£ out ), allowing repetition of site indices in ( T2"5l l and 
bearing in mind that n(£ fl , k/,) + Tl(kb, k a ). 

VII. Discussion 

In this paper we focused on how to generate numerically 
tailored random graphs with controlled macroscopic structural 
properties, to serve e.g. as null models in hypothesis testing. 
Bias in the generation of random graphs has the potential 
to invalidate all further statistical analysis performed on the 
generated networks and has been well documented in the 
literature; it is known to affect the 'stubs' method and the 
'accept-all' edge swapping method. However, the lack so far 
of workable corrections or alternatives has meant that these 
issues have often been ignored. Our theory offers a practical 
and theoretically sound approach to uniformly generating ran- 
dom graphs from ensembles which share certain topological 
characteristics with a real network, and can therefore serve as 
a reliable tool for building unbiased null models. 
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